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The fundamental difference between shear alpha- viscosity and 
turbulent magnetorotational stresses 
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ABSTRACT 

Numerical simulations of turbulent, magnetized, differentially rotating flows driven by the 
magnetorotational instability are often used to calculate the effective values of alpha viscosity 
that is invoked in analytical models of accretion discs. In this paper we use various dynami- 
cal models of turbulent magnetohydrodynamic stresses, as well as numerical simulations of 
shearing boxes, to show that angular momentum transport in MRI-driven accretion discs can- 
not be described by the standard model for shear viscosity. In particular, we demonstrate that 
turbulent magnetorotational stresses are not linearly proportional to the local shear and vanish 
identically for angular velocity profiles that increase outwards. 

Key words: black hole physics - accretion, accretion discs - MHD - instability - turbulence. 



1 INTRODUCTION 

It has long been recognized that molecular viscosity cannot 
be sol ely responsible for angular momentum transport in accretion 
discs. Ishakura & Sunvaevl 1 11 9731) offered an appealing solution to 
this problem by postulating a source of enhanced disc viscosity 
due to turbulence and magnetic fields. The standard accretion disc 
model rests on the idea that the stresses between adjacent disc an- 
nuli are proportional to the local shear, as in a Newtonian laminar 
shear flow, but that it is the interaction of large turbulent eddies that 
results in efficient transport. The idea that turbulent angular mo- 
mentum transfer in accretion discs can be described in terms of an 
enhanced version of the molecular transport operating in (laminar) 
differentially rotating media has been at the core of the majority 
of studies in accretion disc theory an d phenomenology ever since 
(see, e.g. jFrank. King. & Rainell2002l) . 

The origin of the turbulence that leads to enhanced angu- 
lar momentum transport in a ccretion discs has b e en a matter 
of debate since the work of Ishakura & Sunvaevl dl973l) . The 
issue of whether hydrodynamic turbulence can be gener- 
ated and sustained in astrophysical discs, due to the large 
Reynolds number s involved, is currently a matter of re- 
newed interest dAfshordi. Mu khopadhva y. & Naravanl 1 20051 ; 
iMukhopadhvav. Afshordi. & NaravarJ 120051) . However, this 
idea ha s long been chall e nged by analytical dRvu & Goodman! 
1992 1 ; iBalbus & Hawlevl 1200(1, numerical dStone & Balbusl 



lHawlev. Balbus. & Winters!! 19991). and, more rece ntly, experimen- 
tal work ( J i. Burin. Sc hartman. & Goodmarj|2006h . 

During the last decade, it has become evident that the interplay 
between turbulence and magnetic fields is at a more fundamental 
level than originally conceived. There is now strong theoretical 
and numerical evidence suggesting that the process driving the 
turbulence in accretion discs is related to a magnetic instability that 
operates in the presence of a radially decreasing angular velocity 
profile. Since the appreciation of the relevance of this magnetoro- 
tational instability (MRI) to accretion physics (Balbus & Hawley 
1991; see also Balbus & Hawley 1998 and Balbus 2003 for a more 
recen t re view), a variety of local i Hawley. Gammie. & B albus 



I l995l Il996l; Istone. Hawley. Gammie & Balbusl Il99fj: 



[ Brandenburg. Nordlund. Stein. & Torkelssorf 
l 200ll ; ISano. Inutsuka. Turner. & Stonej 



, .mlt . 

2004b 



Brandenburg! 



19961 IBalbus. Hawley. & Stonefll996t IBalbus & Hawlevl 1 19971 ; 



E-mail:mpessah@ias.edu (MEP) 



and global 

dArmitagd Il998l; lHawlevI 120001 . l200ll; lHawlev & Krolikl l200ll; 
IStone & Pringld 120011 numerical simulations have revealed that 
its long-term evolution gives rise to a turbulent state and provides 
a natural avenue for vigorous angular momentum transport. 

The fact that the overall energetic properties of turbulent mag- 
netohydrodynamic (M HD) accretion discs are sim ilar to those of 
viscous accretion discs dBalbus & Papaloiz oul 1 9991) has lead to the 
notion that angular momentum transport due to MRI-driven tur- 
bulence in rotating shearing flows can be described in terms of 
the alpha model proposed by IShakura & Sunvaevl < fl973h . This, in 
turn, has motivated many efforts aimed to comput ing effective al- 
pha values from nu merical simulations (see, e.g.. iGamrni el 1 19981; 
[Brandenburg 1998, and references therein) in order to use them in 
large scale analytical models of accretion discs. 

In this paper, we address in detail how the transport of angular 
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momentum mediated by MHD turbulence depends on the magni- 
tude of the local disc shear and contrast this result to the standard 
model for shear viscosity. We find that one of the fundamental as- 
sumptions in which the standard viscous disc model is based, i.e., 
that angular momentum transport is linearly proportional to the lo- 
cal shear, is not appropriate for describing turbulent MRI-driven 
accretion discs. 



2 ALPHA VISCOSITY VS. TURBULENT MHD 
STRESSES 

The equation describing the dynamical evolution of the mean 
angular momentum density of a fluid element, I, in an axisymmet- 
ric, turbulent MHD accretion disc with tangled magnetic fields is 



dl 1 d , - 

"57 + -7T\ rlvT 
at r or 



r Or 



(1) 



Here, the over-bars denote properly averaged valu es (see, e.g., 
iBalbus & Hawlevll 19981 ; [Balbus & Papaloizoulll999h . v r is the ra- 
dial mean flow velocity, and the quantity T r< j, represents the total 
stress 



T rt h = R rt h — M, 



r4> 



(2) 



acting on a fluid element as a result of the correlated fluctuations in 
the velocity and magnetic fields in the turbulent flow, i.e., 



Rrtf, 



47T 



(3) 
(4) 



where R rt f, and M r( p stand for the r^-components of the Reynolds 
and Maxwell stress tensors, respectively. 

Equation (fTJ highlights the all-important role played by corre- 
lated velocity and magnetic-field fluctuations in turbulent accretion 
discs; if they vanish, the mean angular momentum density of a fluid 
element is conserved. In order for matter in the bulk of the disc to 
accrete, i.e., to loose angular momentum, the sign of the mean to- 
tal stress must be positive. Note that the potential for correlated 
kinetic fluctuations to transport angular momentum in unmagne- 
tized discs is still present in the hydrodynamic version of equation 
(0. However, the dynamical role played by the correlated velocity 
fluctuations in hydrodynamic flows is radically different from the 
corresponding role played by correlated velocity and magnetic field 
fluctuation s in MHD flows, even when th e magnetic fields involved 
are weak (Balb us & Hawievlll997l . ll998h . 

In order to calculate the structure of an accretion disc for 
which angular momentum flows according to equation (T), it is 
necessary to obtain a closed system of equations for the second- 
order correlations defining the total stress T rt f,. In this con- 
text, the original proposal b y Shakura and Sunyaev (see also, 
Lvnden-Bell & Pringle] 1 19741) can be seen as a simple closure 
scheme for the correlations defining the total turbulent stress in 
terms of mean flow variables (e.g., the pressure). 

The model for angular momentum transport on which the 
standard accretion disc is based rests on two distinct assumptions. 
First, it is postulated that the vertically averaged stress exerted on 
any given disc annulus can b e modeled as a shear viscous stress 
dLvnden-Bell & Prinrid 19741) . i.e., in cylindrical coordinates, 



T r4> = -rE!/ turb — . 



(5) 



where E and £1 stand for the vertically integrated disc density and 
the angular velocity at the radius r, respectively. This is a modifica- 



tion of the Newtonian model for the viscous stre ss between adjacent 
layers in a differentially rotating laminar flow dLandau & Lifshita 
1959); in this case, the coefficient fturb parametrizes the turbulent 
kinematic viscosity. In this model, the direction of angular momen- 
tum transport is always opposite to the angular velocity gradient. 
This is the essence of a shear-driven viscous disc. 

Second, on dimensional grounds, it is assumed that the vis- 
cosity coefficient can be parametrized as twb = ac s H. This is 
because the physical mechanism that allows for angular momen- 
tum transport is envisaged as the result of the interaction of tur- 
bulent eddies of typical size equal to the disc scale-height, H, on 
a turnover time of the order of H/c s , where c s is the isothermal 
sound speed. The parameter alpha is often assumed to be constant 
and smaller than unity. In the presence of a shear background, the 
vertically integrated stress is then given by 



rpV 



-aP 



.dlnfi 
dlnr 



(6) 



where P — £cf/7 stands for the average pressure, 7 is the ratio 
of specific heats, and we have used the fact that the scale-height 
of a thin disc in vertical hydrostatic equilibrium is roughly H ~ 
c s /fi. This parametrization of the coefficient of turbulent kinematic 
viscosity implies that the efficiency with which angular momentum 
is transported is proportional to the local average pressure. This is 
the idea behind the standard alpha-disc model. 

It is worth mentioning that the expression usually employed 
to define the stress in alpha-models, i.e., TJL = aP, only provides 
the correct order of magnitude for the stress in terms of the pressure 
for a Keplerian discs. Indeed, in this case, the shear parameter 



dlnQ, 
dlnr 



(7) 



is equal to 3/2 ~ 1. However, the fact that the viscous stress is pro- 
portional to the local shear cannot be overlooked in regions of the 
disc where the local angular velocity can differ significantly from 
its Keplerian valufl This is expected to be the case in the bound- 
ary layer around an accreting star or close to the marginally stable 
orbit around a black hole. These inner disc regions are locally char- 
acterized by different, and possibly negative, shearing parameters 
q. More importantly, if the explicit dependence of the stress on the 
local shear is not considered then equation @ predicts unph ysical, 
non- v anishing stresses for solid body rotation, q = (c.f.. lBlaesl 
120041) . 

More than a decade after the paper by Balbu s~& Hawlevl 

dl99ll) . the MRI stands as the most promising driver of the tur- 
bulence thought to enable the accretion process. The stresses asso- 
ciated with this MRI-driven turbulence have long been considered 
as the physical mechanism behind the enhanced turbulent viscosity 
postulated in the standard model for shear-driven angular momen- 
tum transport. It is important to note, however, that there is no a 
priori reason to assume that the correlated fluctuations defining the 
total turbulent stress in MRI-driven turbulence are linearly propor- 
tional to the local shear. In fact, this assumption can be challenged 
on both theoretical and numerical grounds. 



1 Also note that it is the explicit linear dependence of the stress on the local 
shear what gives the equation for angular momentum transport its diffusive 
nature in standard accretion disc models. 
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3 PREDICTIONS FROM STRESS MODELLING 

There are currently few models that aim to describe the local 
dynamics of turbulent stresses in differentially ro tating magnetized 
media by means of high-ord e r closure schem es (Kato & Yos hizawal 
1 19931 1 19951 : lOgilviel 120031 ; IPessah. Chan. & Psaltisl l2006bH 2 l In 
these models, the total stress, T r< f, = R r< f, — M r $, is not prescribed, 
as in equation {6]l, but its value is calculated considering the local 
energetics of the turbulent flow. This is achieved by deriving a set 
of non-linear coupled equations for the various components of the 
Reynolds and Maxwell stress tensors. These equations involve un- 
known triple-order correlations among fluctuations making neces- 
sary the addition of ad hoc closure relations. 

Although the available models differ on the underlying phys- 
ical mechanisms that drive the turbulence and lead to saturation, 
some important characteristics of the steady flows that they de- 
scribe are qualitatively similar. In particular, all of the models pre- 
dict that turbulent kinetic/magnetic cells in magnetized Keplerian 
discs are elongated along the radial/azimuthal direction, i.e., R rr > 
Rtjytp while M^tp > M rr . Furthermore, turbulent angular momen- 
tum transport is mainly carried by correlated magnetic fluctuations, 
rather than by their kinetic counterpart, i.e., —M r ^ > Rr<t>- All of 
these properties are in general agreement with local numerical sim- 
ulations. However, a distinctive quantitative feature of these models 
that concerns us here is that they make rather different predictions 
for how the stresses depend on the magnitude and sign of the angu- 
lar velocity gradient. 

In the remaining of this section, we highlight the most relevant 
physical properties characterizing the various models and assess 
the functional dependence of the total stress, T rt f,, on the sign and 
steepness of the local angular velocity profile. For convenience, we 
summarize in Appendix A the various sets of equations that define 
each of the models. 



3.1 Kato & Yoshizawa 1995's Model 

In a series of papers, iKato & Yoshizawal d 19931 1 19951) de- 
veloped a model for hydromagnetic tu rbulence in accretion discs 
with no large scale magnetic fields (see lKato. Fukue. & Mineshige 
[1991 for a review, and Nakao 1997 for the inclusion of large 
scale radial and toroidal fields). In their closure scheme, triple- 
order correlations among fluctuations in the velocity and mag- 
netic fields are modeled in terms of second-order correlations using 
the two-scale direct interactio n approximation ( Yos hizawal 1 1 9851 : 
lYoshizawa, Itoh, & Itohll2003l) . as well as mixing l engt h concepts. 

In the model described in lKato & Yoshizawa dl995l) . the shear 
parameter q appears explicitly only in the terms that drive the alge- 
braic growth of the turbulent stresses. The physical mechanism that 
allows the stresses to grow is the shearing of magnetic field lines. 
This can be seen by ignoring all the terms that connect the dynam- 
ical evolution of the Reynolds and Maxwell stresses. When this is 
the case, the dynamics of these two quantities are decoupled. The 
Maxwell stress exhibits algebraic growth while the Reynolds stress 
oscillates if the flow is Rayleigh stable (i.e., q < 2). The growth 
of the magnetic stresses is communicated to the different compo- 



2 The studies described here concern models for the evolution of the var- 
ious correlated fluctuations relevant for describ ing the dynam i cs of a tur- 
bulent magnetized flow. We note that Vishniac & Brandenburg ( 1997) pro- 
posed an incoherent dynamo model for the transport of angular momentum 
driven by the generation of large-scale radial and toroidal magnetic fields. 
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Figure 1. Total dimensionless turbulent stress at saturation as a function of 
the shear parameter q = — dlnQ/dlnr. The various lines show the pre- 
dictions corresponding to the three models for turbulent MHD angular mo- 
mentum transport discussed in [|3] The data points correspond to the volume 
and time averaged MHD stresses calculated from the series of shearing box 
simulations described in ij4] The diamonds and circles represent the results 
obtained in the isothermal (7 = 1) and adiabatic (7 = 5/3) cases, respec- 
tively. The vertical bar in the simulations corresponding to the shear param- 
eter g = 3/2 shows the rms spread in the stresses as calculated from ten 
numerical simulations for a Keplerian disc. This spread, of roughly 20%, 
can be taken as a representative value for the typical rms spread in runs 
with different values of the parameter q for both the isothermal and adi- 
abatic cases. The dotted line corresponds to a linear relationship between 
the stresses and the local shear, like the one assumed in the standard model 
for viscous angular momentum transport. All the quantities in the figure are 
normalized to unity for a Keplerian profile, i.e., for q = 3/2. 

nents of the Reynolds stress via the tensor Sij = Cf R\j — C§ Mij, 
where Cf and Cf are model constants. 

A characteristic feature of this model is that the physical 
mechanism that leads to saturation is conceived as the escape of 
magnetic energy in the vertical direction. This process is incorpo- 
rated phenomenologically by accounting for the leakage of mag- 
netic energy with terms of the form —/3Mij, where [3 stands for the 
escape rate. Although turbulent kinetic and magnetic dissipation 
act as sink terms in the equations for the various stress components 
(proportional to the model parameters ec and em, respectively), if 
it were not for the terms accounting for magnetic energy escape, the 
system of equations would be linear. This means that, in order for a 
stable steady-state solution to exist, either the initial conditions or 
the constants defining the model will have to be fine-tuned. 

The functional dependence of th e total stress on the shear p a- 
rameter q for the model proposed in lKato & Yoshizawal dl995h is 
shown in Figure [Tj The model predicts vanishingly small stresses 
for small values of \q\; the total stress T rr f, is a very strong function 
of the parameter q. Indeed, for q > 0, the model predicts T rt f, ~ q 8 . 



3.2 Ogilvie 2003's Model 

Based on a set of fundamental pr inciples c onstra ining the non- 
linear dynamics of turbulent flows, Ogily3 d2003h developed a 
model for the dynamical evolution of the Maxwell and Reynolds 
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stresses. Five non-linear terms accounting for key physical pro- 
cesses in turbulent media are modeled by considering the form 
of the corresponding triple-order correlations, energy conservation, 
and other relevant symmetries. 

The resulting model describes the development of turbulence 
in hydrodynamic as well as in magnetized non-rotating flows. An 
interesting feature is that, depending on the values of some of the 
model constants, it can develop steady hydrodynamic turbulence 
in rotating shearing flows. For differentially rotating magnetized 
flows with no mean magnetic fields, the physical mechanism that 
allows the stresses to grow algeb raically is the shearin g of magnetic 
field lines, as in the case o f lKato & Yoshizawd (l995). The transfer 
of energy between turbulent kinetic and magnetic field fluctuations 
is mediated by the tensor c 3 M 1/2 M tj - C4.R~ 1/2 MRij, where R 
and M are the traces of the Reynolds and Maxwell stresses, respec- 
tively, and C3 and C4 are model constants. Note that, in this case, the 
terms that lead to communication between the different Reynolds 
and Maxwell stress components are intrinsically non-linear. The 
terms leading to saturation are associated with the turbulent dissipa- 
tion of kinetic and magnetic energy and are given by —ciR 1 ^ 2 Rij 
and —C5M 1 ' 2 Mij, respectively. 

The dependence o f the tota l stress on the shear parameter for 
the model described in lOgilviel J2003I) is shown in Figure Q] For 
angular velocity profiles satisfying < q < 2, the stress be- 
haves as T r< f, ~ q n with n between roughly 2 and 3. Although the 
functional dependence of the total stress on negative values of the 
shear param eter q is different from the one predicted by the model 
described in iKato & Yoshizawal (1995), this model also generates 
MHD turbulence characterized by negative stresses for angular ve- 
locity profiles increasing outwards. 



3.3 Pessah, Chan, & Psaltis 2006's Model 

Motivated by the similarities exhibited by the linear 
regime of the MRI and the fully developed turbulent state 
dPessah. Chan, & Psal tis 2006a), we have recently developed a lo- 
cal model for the growth and saturation of the Reynolds and 
Maxwell stresses in turbulent flows driven by t he magnetorota- 
tional instability dPessah. Chan, & Psaltis! 12006b ) . Using the fact 
that the modes with vertical wave-vectors dominate the fast growth 
driven by the MRI, we obtained a set of equations to describe the 
exponential growth of the different stress components. By propos- 
ing a simple phenomenological model for the triple-order correla- 
tions that lead to the saturated turbulent state, we showed that the 
steady-state limit of the model describes successfully the correla- 
tions a mong stresses found in numerical simulations of shearing 
boxes dHawlev. Gammie. & Balbuslll995h. 

In the model described in IPessah. Chan, & Psalt is (2006b), 
a new set of correlations couple the dynamical evolution of the 
Reynolds and Maxwell stresses and play a key role in develop- 
ing and sustaining the magnetorotational turbulence. In contrast 
to the two previous cases, the tensor connecting the dynamics of 
the Reynolds and Maxwell stresses cannot be written in terms of 
Rij or My. This makes it necessary to incorporate additional dy- 
namical equations for these new correlations. In this model, all the 
second-order correlations exhibit exponential (as opposed to alge- 
braic) growth for shear parameters < q < 2, in agreement with 
numerical simulations. Incidentally, this is the only case in which 
the shear parameter, q, plays an explicit role in connecting the dy- 
namics of the Reynolds and Maxwell stresses. The terms that lead 
to non-linear saturation are proportional to — (M / Mq) 1 / 2 , where 



M is the trace of the Maxwell stress and Mo is a characteristic 
energy density set by the local disc properties. 

The functional dependence of the t otal turbulent stress on 
the loca l shear for the model developed in Pessa h. Chan. & Psaltis 
(2006b) is shown in FigureQ] For angular velocity profiles decreas- 
ing outwards, i.e., for q > 0, the stress behaves like T rt f, ~ q n with 
n between roughly 3 and 4. For angular velocity profiles increas- 
ing outwards, i.e., for q < 0, the stress vanishes identically. Note 
that this is the only model that is characterized by the absence of 
transport for all negative values of the shear parameter q. 

It is worth mentioning that in all three high-order closure 
schemes described above, the pressure, P, does not appear explic- 
itly in the equations defining the models. However, it does play a 
role in setting the overall scale at which the stresses saturate. This 
is because the pressure provides a characteristic velocity (e.g., the 
sound speed) or a characteristic length (e.g., the disc scale height) 
which in turn determine the saturation level of the stresses T r $. In 
order to compare the predictions of how the dimensionless stresses 
depend on the local shear independently of other factors, we nor- 
malized the quantity T r ^/P predicted by each model in Figure Q] 
with the values corresponding to the Keplerian cases with the same 
pressure 0. 

Figure Q] illustrates the sharp contrast between the functional 
dependence of the saturated stresses predicted by all three MHD 
models with respect to the standard shear viscous stress defined in 
equation ((5). It is remarkable that, despite the fact that the various 
models differ on their detailed structure, all of them predict a much 
steeper functional dependence of the stresses on the local shear. 
Indeed, for angular velocity profiles decreasing outwards, they all 
imply T r( p ~ q n with n > 2. The predictions of the various mod- 
els differ more significantly for angular velocit y profiles increasing 
outwa rds. I n this cas e , the m odels developed in lKato & Yoshizawal 
( 1995) and lOgilviel d2003h lead to negative stress es, while the 
model developed in Pessah. Chan. & Psaltis! l l2006bh predicts van- 
ishing stresses. 



4 RESULTS FROM NUMERICAL SIMULATIONS 

There have been only few numerical studies to assess the 
properties of ma gnetorotational turbulence for different valu es of 
the local shear. Abramowicz. Bra ndenburg. & Lasota i ll 9961) car- 
ried out a series of numerical simulations employing the shearing 
box approximation to investigate the dependence of turbulent mag- 
netorotational stresses on the shear-to- vorticity ratio. Although the 
number of angular velocity profiles that they considered was lim- 
ited, their results suggest that the relationship between the turbu- 
lent MHD stresses and the shear is not linear. On the other hand, 
lHawlev. Balbus. & Winters! Jl999n carried out a series of shearing 
box simulations varying the shear parameter from q — 0.1 up to 
q = 1.9 in steps of Ag = 0.1 and reported on the dependence of 



3 Note that the model described in IPessah, Chan. & Psaltij l2006bl) was de- 
veloped to account for the correlations among stresses and magnetic energy 
density at saturation when there is a weak mean magnetic field perpendicu- 
lar to the disc mid-plane. It is known that, for a given initial magnetic energy 
density, numerical simulations of MHD turbulent Keplerian shearing flows 
tend to saturate at higher magn etic energy densities wh en there is a net ver- 
tical magnetic flux lHawlev, Gammie, & Balbus 1996). The normalization 
chosen in Fig. [T] allows us to compare the shear-dependence of the various 
models regardless of their initial field configurations. 
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the Reynolds and Maxwell stresses on the magnitude of the local 
shear but not on the corresponding dependence of the total stress. 

In order to investigate the dependence of MRI-driven turbulent 
stresses on the sign and magnitude of the local shear, we modified 
a version of ZEUS-3D to allow for angular velocity profiles in- 
creasing outwards (i.e., characterized by shear parameters q < 0). 
ZEUS is a publicly available code and is based on an explicit fi- 
nite difference algorithm on a stagger ed mesh. A de tailed descrip- 
tion of this code can be found inlStone & NormarJ jl992al lbh and 
IStone, Mihalas, & Normar] i 19921) . 



4.1 The Shearing Box Approximation 

The shearing box approximation has proven to be fruitful in 
studying the local characteristics of magnetorotational turbulence 
from both the theoretical and numerical points of view. The local 
nature of the MRI allows us to concentrate on scales much smaller 
than the scale height of the accretion disc, H , and regard the back- 
ground flow as essentially homogeneous in the vertical direction. 

In order to obtain the equations describing the dynamics of 
a compressible MHD fluid in the shearing box limit, we consider 
a small box centered at the radius ro and orbiting the central ob- 
ject in corotation with the disc at the local speed vo — ro Qq4>- 
The shearing box approximation consists of a first order expansion 
in r — ro of all the quantities characterizing the flow. The goal of 
this expansion is to retain the most important terms governing the 
dynamics of the MHD flui d in a locally Cartesian coordinate sys- 
tem ( see, e.g.. lGoodman & Xull 19941 : lHawlev. Gammie. & B albus 
This is a good ap proximation as long as t he magnetic fields 
involved are subthermal JPessah & Psa ltis 2005). The resulting set 
of equations is then given by 



carried out up to date (see, e. g., Hawle v. Gammie. & Balbusll 19951: 



- + («.V)« 







(8) 



-V P+ — 

P V 8tt 



(B-V)B 



4irp 

-2to XV + qQlV(r-r ) 2 (9) 

(10) 



^ = V x (v x B) 



dE 
~dt 



+ V ■ (Ev) 



-PV • v 



(11) 



where p is the density, v is the velocity, B is the magnetic field, P 
is the gas pressure, and E is the internal energy density. In writing 
this set of equations, we have neglected the vertical component of 
gravity and defined the local Cartesian differential operator, 



d 



d 



„ _ d 

V = r 7T+ 7n+ z « 
or ro ocp az 



(12) 



where f, <j>, and z are, coordinate-independent, orthonormal ba- 
sis vectors corotating with the background flow at ro. Note that the 
third and fourth terms on the right hand side of equation (|9) account 
for the Coriolis force, present in the rotating frame, and the radial 
component of the tidal field, respectively. We close the set of equa- 
tions <TSt — d 1 It by assuming an ideal gas law with P — (7 — 1)E. 



4.2 Numerical Set Up 

We set the radial, azimuthal, and vertical dimensions of the 
simulation domain to L r = 1, = 6, and L z = 1 and con- 
sider a grid of 32 x 192 x 32 zones. This resolution corresponds 
to the standard resolution used in most shearing box simulations 



ISano. Inutsuka. Turner. & Stonell2004) . 

The density scale in the shearing box is arbitrary and we 
choose po = 1 as in. e.g., [Hawlev. Gammie. & Balbus] d 19951) and 
ISano. Inutsuka. Turner. & Stonel 1 2004 ). We consider the case of 
zero net magnetic flux through the vertical boundariefl by defin- 
ing the initial magnetic field according to B = Bosin[27r(r — 
ro)/L r ]z. The plasma j3 in all the simulations that we perform is 
(3 = P/(Bo/8n) = 200, so the magnetic field is highly subther- 
mal in the initial state. The initial velocity field that corresponds to 
the steady state solution is v = — qflo (r—ro)(p and we choose the 
value do = 10 -3 in order to set the time scale in the shearing box. 
Note that for q = 3/2, this velocity field is simply the first order 
expansion of a steady Keplerian disc around ro. In order to excite 
the MRI, we introduce random perturbations at the 0.1% level in 
every grid point over the background internal energy and velocity 
field in all of the cases. 



4.3 Results 

Keeping all the numerical settings unchanged, we perform 
two suites of numerical simulations with different values of the 
shear parameter q, from q = —1.9 up to q = 1.9 in steps 
of Aq — 0.1. The two sets of runs differ only in the value 
of the adiabatic index 7; we considered an isothermal case, 
with 7 = 1.001, and an adiabatic case, with 7 = 5/3. For 
each value of the shear parameter q, we run each simulation for 
150 orbits. We then compute a statistically meaningful value for 
the saturated stress T r $ and pressure P by averaging the last 
100 orbits in each simulation ((winters. Balbus. & Hawlevl 120031 : 
ISano. Inutsuka. Turner. & Stonell2004h . 

Figure Q] shows the dimensionless stress T r< f,/P obtained for 
both the isothermal and the adiabatic cases (represented with di- 
amonds and circles respectively) as a function of the local shear. 
It is evident from this figure that the turbulent magnetorotational 
stresses are not proportional to the local shear in either the isother- 
mal or the adiabatic cases. There is indeed a strong contrast with 
respect to the standard assumption of a linear relationship between 
the stresses and the local shear (dotted line in the same figure) for 
both positive and negative shear profiles. For angular velocity pro- 
files that decrease with increasing radius, i.e., for q > 0, all of 
the turbulent states are characterized by positive mean stresses and 
thus by outward transport of angular momentum. In these cases, 
stronger shear results in larger saturated stresses but the functional 
dependence of the total stress on the local shear is not linear. For 
angular velocity profiles that increase with increasing radius, i.e., 
for q < 0, all the numerical simulations reach the same final state 
regardless of the magnitude of the shear parameter q. The stresses 
resulting from the initial seed perturbations (at the 0.1% level) 
quickly decay to zero. This is in sharp contrast with the large nega- 
tive stresses that are implied by the standard Newtonian model for 
the shear viscous stress in equations lf5) and (|6). 



4 After submitting this paper we became aware of the impact of numerical 
resolutio n on the saturated stresses in K eplerian sh earing boxes with zero- 
net flu x. IPessah. Ch an. & Psaltis ( 2007j), see also lFromang~& Papaloizou 
120071) . have shown that the saturation of the MRI in this type of simula- 
tions depends linearly on the numerical resolution. The numerical results in 
Figure[T]are normalized with respect to the Keplerian case (q = 1.5). Note 
that the explicit functional dependence of the dimensionless stress on the 
shear, as shown in Figure [T] assumes that resolution affects all the simula- 
tions in the same way, regardless of the value of the shear q. 
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In order to explore further the decay of MHD turbulence 
found for angular velocity profiles increasing outwards, we also 
performed the following numerical experiment. We seeded a run 
with a shear parameter q — —3/2 with perturbations at the 100% 
level of the background internal energy and velocity field. In this 
case, the timescale for the decay of the initial turbulent state was 
longer than the one observed in the corresponding run seeded with 
perturbations at the 0.1% level. The final outcome was nonethe- 
less the same. After a few orbits, the stresses decayed sharply and 
remained vanishingly small until the end of the run at 150 orbits. 
This result highlights the strong stabilizing effects of a positive an- 
gular velocity gradient on the dynamical evolution of the turbulent 
stresses due to tangled magnetic fields. This behavior can be under- 
stood in terms of the joint restoring action due to magnetic tension 
and Coriolis forces acting on fluid elements displaced form their 
initial orbits. 



5 DISCUSSION AND IMPLICATIONS 

In this paper, we investigated the dependence of the turbulent 
stresses responsible for angular momentum transport in differen- 
tially rotating, magnetized media on the local shear as parametrized 
by q = — d In fl/d In r . The motivation behind this effort lies in un- 
derstanding whether one of the fundamental assumptions on which 
much of the standard accretion disc theory rests, i.e., the existence 
of a linear relationship between these two quantities, holds when 
the MHD turbulent state is driven by the MRI. We addressed this 
problem both in the context of current theoretical turbulent stress 
models as well as using the publicly available three-dimensional 
numerical code ZEUS. 

From the theoretical point of view, we have seen that, de- 
spite their differ ent structures, all of the available high-order 
closure schemes dKato & Yoshizawal Il993l 1 19951 ; lOgilviel |2003| ; 
IPessah. Chan. & P saltis 2006bl) predict stresses whose functional 
dependence on the local shear differ significantly from the standard 
model for angular momentum transport. In order to settle this result 
on firmer grounds, we performed a series of numerical simulations 
of MRI-driven turbulence in the shearing box approximation for 
different values of the local shear characterizing the background 
flow. The main conclusion to be drawn from our study is that tur- 
bulent MHD stresses in differentially rotating flows are not linearly 
proportional to the local background disc shear. This finding chal- 
lenges one of the central assumptions in standard accretion disc 
theory, i.e., that the total stress acting on a fluid element in a tur- 
bulent magnetized disc can be modeled as a (Newtonian) viscous 
shear stress. 

We find that there is a strong contrast between the stresses 
produced by MHD turbulence and the viscous shear stresses re- 
gardless of whether the disc angular velocity decreases or increases 
outwards. In the former case, i.e., for q > as in a Keplerian disc, 
the total turbulent stress generated by tangled magnetic fields is not 
linearly proportional to the local shear, q, as it is assumed in stan- 
dard accretion disc theory. On the other hand, for angular velocity 
profiles increasing outwards, i.e., for q < as in the boundary layer 
close to a slowly rotating accreting stellar object, MHD turbulence 
driven by the MRI fails to transport angular momentum, while vis- 
cous shear stresses lead to enhanced negative stresses. 

The functional dependence of the local stress on the shear 
profile determines the topology of transonic accretion flows onto 
black hole s and the radial position of the corresponding criti- 
cal points dAbramowicz & Katolll989l : iKato. Fukue. & Mineshige 



ll998l : lAfshordi & PaczvnskHl2003f) . It also plays a critical role in 
determining the global structure of accretion flows onto stellar ob- 
jects, determining the exchange of angular momentum between the 
disc and the gravitating body and even the angular velocity distri- 
bution itself dPopham & Naravarjfl99lh . If magnetorotational tur- 
bulence is the main mechanism enabling angular momentum trans- 
port in accretion discs then the dependence of the stress on the lo- 
cal shear discussed in this paper can have important implications 
for the global structure and long term evolution of accretion disc 
around proto-stars, proto-neutron stars, accreting binaries, and ac- 
tive galactic nuclei. 
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APPENDIX A: HIGH-ORDER CLOSURE MODELS 

We summarize here the various sets of equations defining the 
models described in Sj3] In order to simplify the comparison be- 
tween the different models, we adopt the notation introduced in Sj2] 
even if this was not the original notation used by the corresponding 
authors. With the same motivation, we work with dimensionless 
sets of equations obtained by using the inverse of the local angular 
frequency (fi^ 1 ) as the unit of time and the relevant characteristics 
speeds or lengths involved in each case. We also provide here the 
values of the various model constants that were adopted in order 
to obtain the total turbulent stresses as a function of the local shear 
shown in Figure Q] 

Al Kato & Yoshizawa 1995' s Model 

The model defined by Kato & Yoshiz awa is based on 
the two-scale direct interactio n approximation (Yosh izawal 19851 : 
lYoshizawa, Itoh. & Itohl 120031) . The temporal evolution of the 
Reynolds and Maxwell stresses is described by 



OtRrr 


— 4R r + n rr — Srr — — £G , 




dtRrtf) 


= (q - 2)R rr + 2R^ + 11,^ - 


S r cj> , 


dtR^tp 


= 2(q — 2)R r £ + Ylq><t> — — 


2 

3 £G 


dtRzz 


= Fizz — 5zz — —EG , 




dtMrr 


2 

= Srr ~ 2(3M r r ~ g £M ' 




d t Mr4> 


= -qM„ + S r4> - 2(3M r4 , , 






= -2qM r4> + S^ - 2f3M u - 


2 

-e M 


d t M zz 


- 2 
= Szz - 2/3M zz - -e M , 
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where the relevant characteristic speed used to define dimension- 
less variables is the local sound speed. 

In this model, the flow of turbulent energy from kinetic to 
magnetic field fluctuations is determined by the tensor 



Si 



C\ Rij 



Co Ma . 



(Al) 



where Cf and C 2 are model constants. This quantity plays the 
most relevant role in connecting the dynamics of the different com- 
ponents of the Reynolds and Maxwell tensors. If Sij is positive, 
the interactions between the turbulent fluid motions and tangled 
magnetic fields enhances the latter. The pressure-strain tensor is 
modeled in the framework of the two scale direct interaction ap- 
proximation according to 

n y = -cF(j^-««iV3)-c?(My-*yM/3) 

— q C^R (SirSjj, + Si^Sjr) , 



(A2) 



where R and M stand for the traces of the Reynolds and Maxwell 
stresses, respectively. This tensor accounts for the redistribution of 
turbulent kinetic energy along the different directions and tends to 
make the turbulence isotropic. The dissipation rates are estimated 
using mixing length arguments and are modeled according to 

ec = ~v G (R + M), (A3) 
3 

e M = ~vm{R + M), (A4) 

where fa and vyi are dimensionless constants. The escape of mag- 
netic energy in the vertical direction is taken into account phe- 
nomenologically via the terms proportional to the (dimensionless) 
rate (3 = XM 1/2 , with < X < 1. 

The values of the constants that we considered in order to ob- 
tain the curv e shown in FigureQJare the same as the ones considered 
in case 2 in lKato & Yoshizawaf l fl995l) . i.e., Co = Cf = Cf = 
Gf = Cl = 0.3, v G = vwi = 0.03. We further consider X = 0.5 
as a representative case. 



A2 Ogilvie 2003's Model 

Ogilvie proposed that the functional form of the equations 
governing the local, non-linear dynamics of a turbulent MHD flow 
are strongly constrained by a set of fundamental principles. The 
model that he developed is given by the following set of equations 

d t R rr = m r4> - Cl R 1/2 R rr - C 2 R 1,2 {Rrr- R/$) 
+C 3 M 1/2 M rr ~ C4,R~ 1/2 MRrr , 

dtRrj, = (q- 2)R rr + 2R 4>4> ~ (ci + c 2 )R 1/2 R r<t> 



+c 3 M 1/2 M r4 , ~ c 4 R- 1/2 MR r4 , , 
dtR^ = 2{q-2)R r < l ,-c 1 R 1/2 R^~c 2 R 1/2 (R^- R/3) 
+c- A M 1/2 M 4>4 , - c 4 R~ 1/2 MR^ , 



d t R z 



= - Cl R 1/2 R 



C2 



R 1/2 (R zz -R/3) 



+c 3 M 1/2 M zz - Ci R- 1/2 MR zz , 
c 4 R- 1/2 MR r 



d t M rr = c 4 R- 1/2 MR rr -{c 3 + cs)M 1/2 M rr , 

d t M r4> = -qM rr + c A Rr 1/2 MR r ^ - (c 3 + c B )M 1/2 M r ^ , 

dtM^ = -2qM r4 , + c 4 R- 1/2 MR 4 , 4 , - (c 3 + c 5 )M 1/2 M M , 

dtM zz = CiR- 1/2 MR zz - (c3 + c 5 )M 1/2 M zz . 



Here R and M denote the traces of the Reynolds and Maxwell 
tensors and we have defined the quantities ci , . . . , cs which are 
related to the positive dimensionless constants defined by Ogilvie 



Ci, . . . , Cg via d = Ld, where L is a vertical characteristic 
length (e.g., the thickness of the disc). Note that Ogilvie's origi- 
nal equations are written in terms of Oort's first constant A — q/2 
(in dimensionless units). 

In this model, the constant c 2 dictates the return to isotropy 
expected to be exhibited by freely decaying hydrodynamic turbu- 
lence. The terms proportional to C3 and a transfer energy between 
kinetic and magnetic turbulent fields. The constants ci and C5 are 
related to the dissipation of turbulent kinetic and magnetic energy, 
respectively. Note that, in order to obtain the representative behav- 
ior of the total turbulent stress as a function of the local shear that 
is shown Figure [T] we set c\ , . . . , C5 = 1. 



A3 Pessah, Chan, & Psaltis 2006's Model 

We have recently developed a local model for the growth and 
saturation of the Reynolds and Maxwell stresses in turbulent flows 
driven by the magnetorotational instability that leads to exponen- 
tial growth for the stresses and can account for a number of corre- 
lations observed in numerical simulations ( Pess ah. Chan. & Psaltisl 
2006b). In this model, the Reynolds and Maxwell stresses are not 
only coupled by the same linear terms that drive the turbulent state 
in the previous two models but there is also a new tensorial quan- 
tity that couples their dynamics further. This new tensor cannot be 
written in terms of Rij or Mij , making it necessary to incorporate 
additional dynamical equations. 

The set of equations defining this model is 



OtRr 



AR r<t , + 2W r<t 




dtRcj, = (q- 2)Rr r + 2R^ - Wrr + 




dtRd 



= 2(q-2)R r4 ,-2W^ r -\ 7r R^, 

Mo 



dtWrr = qW r4> + 2WV + C X„(fir* - M r4> ) - 



22 - - M - 

dtWrj, = 2W^ - C fc max (i?rr ~ M rr ) - \ -^-Wrj, , 

V Mo 




OtW^r = (q-2)W rr +qW U 




M H7 



dtW^ = {q - 2)W r< f> - C A&axCRr* - Mri) 



d t M rr = -2Wrd, - \ \ ITT M r r ■ 




d t M r 



-qMrr + W r 



-2qMr4, + 2W^, r - 4/ -^M^ . 

Mo 




where we have defined dimensionless variables considering the 
mean Alfven speed v\ z = B z /y /r 4npo, with B z the local mean 
magnetic field in the vertical direction and po the local disc den- 
sity. The tensor Wit = (SviSjk) is defined in terms of correlated 
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fluctuations in the velocity and current (Sj = V X SB) fields. The 
(dimensionless) wavenumber defined as 

^max = Q (A5) 

corresponds to the scale at which the MRI-driven fluctuations ex- 
hibit their maximum growth and 

M = £p Hn v A z , (A6) 

is a characteristic energy density set by the local disc properties, 
with H the disc thickness. The parameters £ ~ 0.3 and £ ~ 11 
are model constants which are determined by requiring that the 
Reynolds and Maxwell stresses satisfy the correlations observed in 
numerical simulations of Kepleri an shearing boxes with q — 3/2 
( THawlev. Gammie. & Balbusll995h . 
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